function [output subjectinfo] = TryClustering_AllSubjects(subjectinfo, cfg)

subjects = {'Subject DS', 'Subject CW','Subject DJ','Subject RV'};
%'Subject FG','Subject JF', geen FaceHouse

for s = 1:length(subjects)
    subject=subjects{s};
    clearvars subjectinfo;
    
    if ~exist('subjectinfo', 'var'); subjectinfo = SubjectInfo.Load(subject); end;
    if ~exist('cfg','var'); cfg = TryClustering_DefaultCfg(); end;
    load(subjectinfo.EpochData.(cfg.EpochSet).TfcDatafile);
    
    %% Ff correctie voor tfc.times
    tfc.times = tfc.times - 5000;
    
    datatable = BuildDatatable(tfc, subjectinfo);
    X = cell2mat(datatable(2:end, 5));
    Y = cell2mat(datatable(2:end, 6));
    Z = cell2mat(datatable(2:end, 4));
    
    regions = unique(datatable(2:end, 1));
    bins = unique(Z);
    clusters = [];
    for r = regions'
        for b = bins'
            bintf = Z == 1;
            regtf = strcmp(cell2mat(r), datatable(2:end, 1));
            dataValues = Y(bintf & regtf);
            dataTimes = tfc.times(X(bintf & regtf));
            clusters = [clusters FindClusters(dataTimes, dataValues, tfc.subject, tfc.epochs, b, cell2mat(r),cfg)];
        end
    end
    
    % Teken het spul..
    
    % figure; hold on; grid on;
    % for c = 1:length(clusters)
    %     rx = max(abs([clusters(c).members.time]-clusters(c).meanTime))*1.5;
    %     ry = max(abs([clusters(c).members.value]-clusters(c).meanValue))*1.5;
    %     circle(clusters(c).meanTime, clusters(c).meanValue, [rx ry]); %, [0.95 0.95 0.95], [1 0.5 0.5], 0, 360, 0, 2);
    %     scatter([clusters(c).members.time], [clusters(c).members.value], 40, 'filled');
    %     % ezcontour(@(x,y)pdf(obj,[x y]),[min(times(members))-200 max(times(members))+200], [min(vals(members))-0.2 max(vals(members))+0.2]);
    % end
    
    % output.clusters = clusters;
    
end
end


function datatable = BuildDatatable(tfc, subjectinfo)
datatable = {'Region', 'Electrode', 'Bin', 'BinNr', 'Latency', 'Value'};
electrodes = fieldnames(tfc.data);
for i=1:length(electrodes)
    electrode = cell2mat(electrodes(i));
    currentRow.Electrode = electrode;
    currentRow.Region = ElectrodeInfo.RegionByLabel(subjectinfo.Electrodes, electrode);
    for j=1:length(tfc.data.(electrode));
        d = tfc.data.(electrode)(j);
        % Find the bin description
        currentRow.BinNr = find(tfc.freqbins(:,1)==d.frequencyrange(1) & tfc.freqbins(:,2)==d.frequencyrange(2));
        currentRow.Bin = tfc.freqbindescriptions{currentRow.BinNr};
        for p = 1:length(d.peaklocs)
            currentRow.Latency = d.peaklocs(p);
            currentRow.Value = d.mean(d.peaklocs(p));
            % Write the data to the table..
            datatable(end+1,:) = { ...
                currentRow.Region, ...
                currentRow.Electrode, ...
                currentRow.Bin, ...
                currentRow.BinNr, ...
                currentRow.Latency, ...
                currentRow.Value};
        end
    end
end
end

function [clusters] = FindClusters(times, vals, subject, epochset, bin, region ,cfg)
clusters = [];
[T I] = sort(times)
times = times(I)
vals = vals(I)
membership = ones(length(times),1)*NaN;
nextCluster = 1;
for i=1:length(times)
    if isnan(membership(i));
        membership(i) = nextCluster;
        nextCluster = nextCluster+1;
    end
    for j=i:length(times)
        % If the distance in time is <300 and the type of peak (Neg or Pos)
        % is the same, they belong to the same cluster.
        times(i)-times(j)
        abs(times(i)-times(j))
        if abs(times(i) - times(j))<300 && sign(vals(i)-1) == sign(vals(j)-1);
            membership(j) = membership(i);
        end
    end
end

for c = min(membership):max(membership)
    members = find(membership==c);
    
    if length(members)<2; continue; end;
    
    tcluster= max(times(members))-min(times(members));
    lambda=length(times)/(8000/tcluster);
    kcluster=length(members);
    
    chance = (exp(-lambda)*lambda^kcluster)/factorial(kcluster);
    
    if chance>0.05; continue; end;
    
    %     if weight<0.8; continue; end;
    
    %      weight = length(members)/electr; %% ratio of # of point in cluster and # of electrodes in the region
    
    Tm = mean(times(members));
    Vm = mean(vals(members));
    Vstd = std(vals(members));
    Tstd = std(times(members));
    
    %    scatter(times(members),vals(members),40,'filled');
    
    mCfgMs = fieldnames(cfg.Markers);
    markerIdx=[];
    
    for i = 1:length(mCfgMs)
        if cfg.Markers.(mCfgMs{i}).Description == region
            markerIdx = mCfgMs{i};
            break;
        end
        if isempty(markerIdx); markerIdx = 'Other'; end;
    end
    
    markerType = cfg.Markers.(markerIdx).Marker;
    markerEdgeSize = cfg.Markers.(markerIdx).MarkerEdgeSize;
    markerSize = cfg.Markers.(markerIdx).MarkerSize;
    markerEdgeColor = cfg.Markers.(markerIdx).MarkerEdgeColor;
    markerFaceColor = cfg.Markers.(markerIdx).MarkerFaceColor;
    
    herrorbar(Tm,Vm,Tstd);
    errorbar(Tm,Vm,Vstd,...
        'Marker', markerType, ...
        'MarkerSize', markerSize, ...
        'MarkerEdgeColor', markerEdgeColor, ...
        'LineWidth',0,...
        'MarkerFaceColor', markerFaceColor); hold on;
    grid on;
    
    clearvars markerIdx;
    
    currentCluster = c - min(membership) + 1; % start from 1
    
    clusters(currentCluster).meanTime = mean(times(members));
    clusters(currentCluster).meanValue = mean(vals(members));
    clusters(currentCluster).stdTime = std(times(members));
    clusters(currentCluster).stdValue = std(vals(members));
    
    for i = members'
        clusters(currentCluster).members(i).time = times(i);
        clusters(currentCluster).members(i).value = vals(i);
    end
    % add distribution
    %    D = [times(members)' vals(members)];
    %    [n, d] = size(D);
    %   if n > d
    %       clusters(currentCluster).gmdistribution = gmdistribution.fit(D,1);
    %   else
    %       clusters(currentCluster).gmdistribution = [];
    %    end
    % add other info
    
    clusters(currentCluster).subject = subject;
    clusters(currentCluster).epochset = epochset;
    clusters(currentCluster).bin = bin;
    clusters(currentCluster).region = region;
    
end
end
